\defmodule {NormalInverseGaussianProcess}

This class represents a normal inverse gaussian process (NIG).
It obeys the stochastic differential equation \cite{pBAR98a}
\begin{equation}
   dX(t) = \mu dt + dB(h(t)), \label{eq:nig}
\end{equation}
where $\{B(t),\, t\ge 0\}$ is a \class{BrownianMotion} 
with drift $\beta$ and variance 1,
and $h(t)$ is an \class{InverseGaussianProcess} $IG(\nu/\gamma,\nu^2)$, with 
$\nu = \delta dt$ and $\gamma = \sqrt{\alpha^2 - \beta^2}$.

In this class, the process is generated using the sequential
technique:  $X(0)=x_0$ and 
\begin{equation}
   X(t_j) - X(t_{j-1}) =\mu dt +  \beta Y_j + \sqrt{Y_j} Z_j, 
\end{equation}
where $Z_j \sim N(0,1)$, and $Y_j \sim IG(\nu/\gamma,\nu^2)$
with $\nu = \delta(t_j - t_{j-1})$.   

There is one \externalclass{umontreal.iro.lecuyer.rng}{RandomStream} 
used to generate the $Z_j$'s and
there are one or two streams used to generate the underlying
\class{InverseGaussianProcess}, depending on which IG subclass
is used.

In finance, a NIG process usually means that
the log-return is given by a NIG process;
\class{GeometricNormalInverseGaussianProcess}
should be used in that case.


\bigskip\hrule\bigskip

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{code}
\begin{hide}
/*
 * Class:        NormalInverseGaussianProcess
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
\end{hide}
package umontreal.iro.lecuyer.stochprocess;\begin{hide}
import umontreal.iro.lecuyer.rng.*;
import umontreal.iro.lecuyer.probdist.*;
import umontreal.iro.lecuyer.randvar.*;
import umontreal.iro.lecuyer.hups.*;

\end{hide}

public class NormalInverseGaussianProcess extends StochasticProcess \begin{hide} {
    protected RandomStream streamIG1;   // U[0,1) gen used in the inverse gaussian
    protected RandomStream streamIG2;   // U[0,1) gen used (maybe) in the inverse gaussian
    protected RandomStream streamBrownian;   // U[0,1) gen for the normal "Brownian"

    protected InverseGaussianProcess igProcess; 
    protected NormalGen normalGen;
    // Randomized time generated by the InverseGaussianProcess.
    protected double[] stochTime;
    double[] dt;
    double[] mudt;

    // Parameters of the normal inverse gaussian
    protected double mu;
    protected double delta;
    protected double alpha;
    protected double beta;
    protected double gamma;


\end{hide}
\end{code}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Constructors}
\begin{code}

   public NormalInverseGaussianProcess (double x0, double alpha,
                                        double beta, double mu,
                                        double delta,
                                        RandomStream streamBrownian,
                                        InverseGaussianProcess igP) \begin{hide} {
        setParams (x0, alpha, beta, mu, delta);
        this.streamBrownian = streamBrownian;
        normalGen = new NormalGen(streamBrownian); // N(0,1)
        igProcess = igP;  // its params will be set in init().
        this.streamIG1 = igProcess.getStream();
        this.streamIG2 = streamIG1;
    }\end{hide}
\end{code}
\begin{tabb} Given an \class{InverseGaussianProcess} \texttt{igP}, constructs a
new \texttt{NormalInverseGaussianProcess}.  The parameters and observation
times of the IG process will be overriden by the parameters
of the NIG process.  If there are two 
\externalclass{umontreal.iro.lecuyer.rng}{RandomStream}'s in the 
\class{InverseGaussianProcess}, this constructor assumes that
both streams have been set to the same stream.
\end{tabb}
\begin{code}

   public NormalInverseGaussianProcess (double x0, double alpha,
                                        double beta, double mu,
                                        double delta,
                                        RandomStream streamBrownian,
                                        RandomStream streamIG1,
                                        RandomStream streamIG2,
                                        String igType) \begin{hide} {
        setParams (x0, alpha, beta, mu, delta);
        this.streamIG1 = streamIG1;
        this.streamIG2 = streamIG2;
        this.streamBrownian = streamBrownian;
        normalGen = new NormalGen(streamBrownian); // N(0,1)

        // The initial time of igProcess here, 0., is arbitrary: 
        // init() sets igProcess.x0 = t[0].
        if(igType.compareTo("SEQUENTIAL_SLOW") == 0)
            // the initial value, 0.0, of the IG will be overriden by
            // the initial time, when it will be set.
            igProcess = new InverseGaussianProcess(0.0, delta, gamma, 
                                streamIG1);
        else if(igType.compareTo("SEQUENTIAL_MSH") == 0)
            igProcess = new InverseGaussianProcessMSH (0.0, delta, gamma, 
                                streamIG1, streamIG2);
        else if(igType.compareTo("BRIDGE") == 0)
            igProcess = new InverseGaussianProcessBridge (0.0, delta, gamma, 
                                streamIG1, streamIG2);
        else if(igType.compareTo("PCA") == 0)
            igProcess = new InverseGaussianProcessPCA (0.0, delta, gamma, 
                                streamIG1);
        else throw new IllegalArgumentException("Unrecognized igType");
    }\end{hide}
\end{code}
\begin{tabb} 
Constructs a new \texttt{Normal\-Inverse\-Gaussian\-Process}.
The string argument corresponds to the type of underlying
\class{Inverse\-Gaussian\-Process}.  The choices are SEQUENTIAL\_SLOW,
SEQUENTIAL\_MSH, BRIDGE and PCA, which correspond
respectively to \class{Inverse\-Gaussian\-Process}, \class{Inverse\-Gaussian\-ProcessMSH},
\class{Inverse\-Gaussian\-Process\-Bridge} and \class{Inverse\-Gaussian\-ProcessPCA}.
The third \externalclass{umontreal.iro.lecuyer.rng}{Random\-Stream}, \texttt{streamIG2},
will not be used at all if the SEQUENTIAL\_SLOW or PCA methods are chosen.
\end{tabb}
\begin{code}

   public NormalInverseGaussianProcess (double x0, double alpha,
                                        double beta, double mu,
                                        double delta,
                                        RandomStream streamAll,
                                        String igType) \begin{hide} {
        this(x0, alpha, beta, mu, delta, streamAll, streamAll, streamAll, igType);
    }\end{hide}
\end{code}
\begin{tabb} Same as above, but all \externalclass{umontreal.iro.lecuyer.rng}{RandomStream}'s 
are set to the same stream, \texttt{streamAll}.
\end{tabb}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Methods}
\begin{code}

   public double[] generatePath() \begin{hide} {
        if(igProcess.getNumberOfRandomStreams() != 1)
            return generatePathTwoIGStreams();

        double x = x0;
        double[] randomNormal = new double[d];
        double[] randomIG1 = new double[d];
        for (int j = 0; j < d; j++)
        {
            randomIG1[j]    = streamIG1.nextDouble();
            randomNormal[j] = streamBrownian.nextDouble();
        }

        stochTime = igProcess.generatePath(randomIG1);
        for (int j = 0; j < d; j++)
        {
            double dy = stochTime[j + 1] - stochTime[j];
            x += mudt[j] + beta*dy + 
                 Math.sqrt(dy) * NormalDistQuick.inverseF01(randomNormal[j]);
            path[j + 1] = x;
        }
        observationIndex = d;
        return path;
    }\end{hide}
\end{code}
\begin{tabb} 
Generates the path.  
This method samples each stream alternatively,
which is useful for quasi-Monte Carlo, where
all streams are in fact the same iterator on 
 a \externalclass{umontreal.iro.lecuyer.hups}{PointSet}.
\end{tabb}
\begin{code}\begin{hide}

    protected double[] generatePathTwoIGStreams() {
        double x = x0;
        double[] uniformNormal = new double[d];
        double[] uniformIG1 = new double[d];
        double[] uniformIG2 = new double[d];

        for (int j = 0; j < d; j++)
        {
            uniformIG1[j]    = streamIG1.nextDouble();
            uniformNormal[j] = streamBrownian.nextDouble();
            uniformIG2[j]    = streamIG2.nextDouble();
        }

        stochTime = igProcess.generatePath(uniformIG1, uniformIG2);
        for (int j = 0; j < d; j++)
        {
            double dy = stochTime[j + 1] - stochTime[j];
            x += mudt[j] + beta*dy + 
                 Math.sqrt(dy) * NormalDistQuick.inverseF01(uniformNormal[j]);
            path[j + 1] = x;
        }
        observationIndex = d;
        return path;
    }\end{hide}

   public double nextObservation() \begin{hide} {
       double igNext = igProcess.nextObservation();
       observationIndex = igProcess.getCurrentObservationIndex();
       stochTime[observationIndex]  = igNext;
       double dY = igNext - stochTime[0];
       double dT = t[observationIndex] - t[0];
       path[observationIndex] = x0  +  mu * dT  +  beta * dY  +
                                Math.sqrt(dY) * normalGen.nextDouble();
       return path[observationIndex];
    }\end{hide}
\end{code}
\begin{tabb} 
Returns the value of the process for the next time step.
If the underlying \class{InverseGaussian\-Process}
is of type \class{InverseGaussianProcessPCA}, this method cannot be
used.  It will work with \class{InverseGaussianProcessBridge}, but
the return order of the observations is the bridge order.
\end{tabb}
\begin{code}\begin{hide}

    protected void init()
    {
        super.init();
        igProcess.setParams(delta, gamma);
        if(observationTimesSet){
            stochTime = new double[d+1];
            stochTime[0] = t[0];
            dt = new double[d];
            mudt = new double[d];
            for(int i = 0; i < d; i++){
               dt[i]   = t[i+1] - t[i];
               mudt[i] = dt[i] * mu;
            }
        }
    }\end{hide}

   public void setObservationTimes(double t[], int d) \begin{hide} {
        super.setObservationTimes(t, d);
        igProcess.setObservationTimes(t, d);
        igProcess.x0 = t[0];
    }\end{hide}
\end{code}
\begin{tabb} Sets the observation times on the NIG process as usual,
but also sets the observation times of the underlying \class{InverseGaussianProcess}.
It furthermore sets the starting \emph{value} of the \class{InverseGaussianProcess}
to \texttt{t[0]}.
\end{tabb}
\begin{code}

   public void setParams (double x0, double alpha, double beta,
                          double mu, double delta) \begin{hide} {
        // Initializes the parameters of the process
        if (delta <= 0.0)
            throw new IllegalArgumentException ("delta <= 0");
        if (alpha <= 0.0)
            throw new IllegalArgumentException ("alpha <= 0");
        if (Math.abs(beta) >= alpha)
            throw new IllegalArgumentException ("|beta| >= alpha");

        this.x0    = x0;
        this.mu    = mu;
        this.delta = delta;
        this.beta  = beta;
        this.alpha = alpha;

        gamma = Math.sqrt(alpha*alpha - beta*beta);
        if (observationTimesSet) init();
    }\end{hide}
\end{code}
\begin{tabb} Sets the parameters. Also, computes 
$\gamma = \sqrt{\alpha^2-\beta^2}$.
\end{tabb}
\begin{code}

   public double getAlpha() \begin{hide} {
      return alpha;
   }\end{hide}
\end{code}
\begin{tabb} Returns alpha.
\end{tabb}
\begin{code}

   public double getBeta() \begin{hide} {
      return beta;
   }\end{hide}
\end{code}
\begin{tabb} Returns beta.
\end{tabb}
\begin{code}

   public double getMu() \begin{hide} {
      return mu;
   }\end{hide}
\end{code}
\begin{tabb} Returns mu.
\end{tabb}
\begin{code}

   public double getDelta() \begin{hide} {
      return delta;
   }\end{hide}
\end{code}
\begin{tabb} Returns delta.
\end{tabb}
\begin{code}

   public double getGamma() \begin{hide} {
      return gamma;
   }\end{hide}
\end{code}
\begin{tabb} Returns gamma.
\end{tabb}
\begin{code}

   public double getAnalyticAverage (double time) \begin{hide} {
        return mu*time+delta*time*beta/gamma;
    }\end{hide}
\end{code}
\begin{tabb} Returns the analytic average, which
is $\mu t + \delta t \beta/ \gamma$.
\end{tabb}
\begin{code}

   public double getAnalyticVariance (double time) \begin{hide} {
        return delta*time*alpha*alpha/gamma/gamma/gamma;
    }\end{hide}
\end{code}
\begin{tabb} Returns the analytic variance, which is
$\delta  t  \alpha^2 / \gamma^3$.
\end{tabb}
\begin{code}

   public RandomStream getStream() \begin{hide} {
        if( (streamIG1 != streamIG2)                ||
            (streamIG1 != streamBrownian)         ||
            (streamIG1 != normalGen.getStream())  ||
            (streamIG1 != igProcess.getStream())  )
            throw new UnsupportedOperationException
            ("Two different streams or more are present");
        return streamIG1;
    }\end{hide}
\end{code}
\begin{tabb} Only returns the stream if all streams are equal,
including the stream(s) in the underlying \class{InverseGaussianProcess}.
\end{tabb}
\begin{code}

   public void setStream (RandomStream stream) \begin{hide} {
        streamIG1 = streamIG2 = streamBrownian = stream;
        normalGen.setStream(stream);
        igProcess.setStream(stream);
    }\end{hide}
\end{code}
\begin{tabb} Sets all internal streams to \texttt{stream},
including the stream(s) of the underlying \class{Inverse\-Gaussian\-Process}.
\end{tabb}
\begin{code}
\begin{hide}
} 
\end{hide}
\end{code}
